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Abstract. In this work, we present two numerical methods to approximate 
solutions of systems of dissipative sine-Gordon equations that arise in the study 
of one-dimensional, semi-infinite arrays of Josephson junctions coupled through 
superconducting wires. Also, we present schemes for the total energy of such 
systems in association with the finite-difference schemes used to approximate 
the solutions. The proposed methods are conditionally stable techniques that 
yield consistent approximations not only in the domains of the solution and 
the total energy, but also in the approximation to the rate of change of energy 
with respect to time. The methods are employed in the estimation of the 
threshold at which nonlinear supratransmission takes place, in the presence 
of parameters such as internal and external damping, generalized mass, and 
generalized Josephson current. Our results are qualitatively in agreement with 
the corresponding problem in mechanical chains of coupled oscillators, under 
the presence of the same parameters. 



1. Introduction 

The well-known sine-Gordon equation is an important differential equation in 
physics due to its multiple applications. For example, it appears in the study of long 
Josephson junctions between superconductors [1], in the investigation of discrete 
arrays of Josephson junctions coupled through superconducting wires [2], in the 
modeling of the motion of a damped string in a non-Hookean medium or, more 
generally, in the description of the motion of rigid pendula attached to a stretched 
wire [3] , and as a model for rapidly rotating baroclinic fluids . Modified versions 
of the sine-Gordon equation often appear in the form of Klein-Gordon equations, 
in problems such as the study of fluxons in Josephson transmission lines [5] , in the 
statistical mechanics of nonlinear coherent structures such as solitary waves (see [5] 
pp. 298-309), in the study of Alfven waves in nonuniform media [7], or as Taylor 
approximations to problems that involve the sine-Gordon model. 

The continuous form of the sine-Gordon equation possesses many interesting ana- 
lytical properties. For example, the sine-Gordon equation is an integrable equation, 
a characteristic which does not share, for instance, with the nonlinear Klein-Gordon 
equation. The existence of traveling-wave solutions, soliton solutions and, more- 
over, nonlinear intrinsic modes (called moving breathers) , are also well-known facts 
associated to this equation [5] . In summary, many interesting properties are already 
available for the analysis of continuous sine-Gordon systems; however, the theory 



2010 Mathematics Subject Classification. (PACS) 45.10.-b; 05.45.-a; 02.30.Hq. 
Key words and phrases, finite-difference schemes; sine-Gordon equation; discrete Josephson- 
junction arrays; nonlinear supratransmission; stability analysis. 



2 



J. E. MACIAS-DIAZ AND S. JEREZ-GALIANO 



underneath the foundations of this differential equation is rather far from being 
completely understood, whence it follows that the determination of novel results 
related to continuous sine-Gordon media, as well as the study of new physical ap- 
plications, are highly transited highways in the literature of applied mathematics 
and physics. 

From a numerical perspective, the study of the classical, continuous sine-Gordon 
equation has been a popular source of interest since its inception [9 . For instance, 
Josephson tunnel junctions were studied numerically via the sine-Gordon equation 
in [S] , while the numerical study of wave transmission in long Josephson structures 
subject to harmonic driving was begun by Olsen and Samuelsen ^10 . On the other 
hand, the (l-l-l)-dimensional Langevin equation was solved through numerical sim- 
ulations in 11 . Later on, Strauss and Vazquez designed a numerical technique to 
approximate radially symmetric solutions of conservative Klein-Gordon equations 
|12) , one of the most important features of their numerical method being that the 
discrete energy associated with the differential equation is conserved. 

It must be mentioned that, after Strauss and Vazquez's pioneering work, the 
development of numerical techniques with properties of invariance in some physical 
domains became a very fruitful topic of research. Particularly, many finite-difference 
schemes (implicit and explicit) with properties of invariance in the energy domain 
were proposed for classes of nonlinear, conservative Klein-Gordon equations |12l 
[T3l [Ml [15] . However, it is worth mentioning that the nonconservative case has 
been left aside, and that its study has been recently initiated [El [171 Ull [19] - 
motivated by the study of the nonlinear phenomenon of supratransmission of energy 
in certain dissipative systems, where the sine-Gordon equation clearly dominates 
the theoretical [H [^Hl [H] and the applied [HI [^^ scenarios. 

This article presents two conditionally stable, finite-difference schemes with prop- 
erties of consistency in the domains of the solution and the energy. The numerical 
techniques are used in the investigation of the process of nonlinear supratransmis- 
sion in discrete arrays of Josephson junctions subject to harmonic perturbations in 
one end, for which the methods presented here seem to be natural choices in the 
study of the problem. 

Section [2] of this paper introduces the mathematical model under study and 
the energy expressions employed in the analysis of supratransmission. Section [3] 
presents two numerical schemes associated to our problem and some practical ob- 
servations on the implementation of the methods. The next section establishes 
the most important numerical properties of the computational techniques; here, we 
show that the proposed methods provide consistent approximations in the energy 
domain — a highly desirable characteristic in the analysis of the process of supra- 
transmission — , and establish stability properties. In Section [5l our techniques are 
applied to the determination of the supratransmission threshold in the presence of 
several parameters. Finally, the article closes with a section of concluding remarks 
and directions of further research. 

2. Preliminaries 

2.1. Mathematical model. The physical motivation of the present work is the 
study of discrete, linear arrays of parallel Josephson junctions coupled through su- 
perconducting wires, subject to harmonic driving on one end. The model considers 
the effects of internal and external damping, relativistic mass, and a generalized 
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Josephson current. More concretely, we consider a sequence of iV coupled Josephson 
junctions initially at rest, and study the associated initial-value problem 

iii - (c^ + ;3A) A^wi -f- 7-U1 + m^ui + sinui = J + 0(i), 

(1) 



Un — (c^ + PDt) A^Un + jiin + rn^7i„ + sinu„ = ,/, 1 < n < N, 
UN + (c^ + POt) A^un-i + JUn + rn^UN + sinujM — J — I, 



( M„(0)=0, l<n<iV, 

s. t. < du„ , , 

I -^(0) =0, \<n<N. 

Here, Dt is the derivative operator with respect to time i, 

(2) A^U„ = (m„-|-1 - Un)X{\,N\ {n + 1) + (Wn-l - Un)xilM] ~ 1) , 

and is the characteristic function on the set [1,-/V] evaluated at n, which 

is equal to 1 if n belongs to [1, A^], and equal to otherwise. 

Notationally, c is a nonnegative constant called the coupling coefficient, and (3 
and 7 are also nonnegative constants called the internal and external damping 
coefficients, respectively. The constant m — called the mass of the system — is 
a real or pure-imaginary complex number that has been included in the system of 
equations in order to suggest possible, further applications of our technique to phi- 
four theories and super symmetry |24j . Meantime, the function is called the input 
intensity function, and it is a harmonic disturbance irradiating at a frequency in the 
forbidden band-gap of the system, that takes the concrete form (j>{t) = Asm{Vtt). 
Finally, the constants J and / are called, respectively, the generalized Josephson 
current and the output current intensity of the system. 

It has been established [25 that the inclusion of two convenient nodes, one 
located at the beginning and the other at the end of the array of the N Josephson 
junctions in ([1]), transforms our initial- value problem into the equivalent initial- 
value problem with boundary data 

Un - (c^ + PDt) Alun + jnUn + Tn^u„ + siu u„ = J, I <n < N, 
( un{0) = 0, M„(0) = 0, l<n<N, 
^' s.t. I c^{uo-ui) + l3{uo-ui) = (j), t>0, 

{ UN+i - UN = 0, t>Q, 

where 7n = 7 for every n < N, and jn = 7 + ^/R- The number R is called 
the output reading resistance of the system, and it is related to the output current 
intensity through Ohm's law: I — un / R- 

Let us consider the case /3 — 7 = J — 0. Then, the linear dispersion relation 
associated with the linearized form of the differential equations in (jS]) is given by 
oP{k) = + 1 + 2c^(l — cos k). It is clear that 1 — cos fc > for every k, so that 
w^(fc) > + It follows that the forbidden band-gap is provided by the inequality 
n < + 1. 

2.2. Energy expressions. It is worth noticing that the Hamiltonian of the n-site 
in the conservative version of ([3]) is given by 

(4) Hn = ^[ul + C^{Un+l - Unf + m^U^] + V (Un) ~ Ju„ , 

with V{un) being the potential function for the classical sine-Gordon equation 
evaluated at u„, namely, V{un) — 1 — cos u„. The inclusion of the potential between 
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the coupling of the first two sites in the chain of Josephson junctions results in the 
following expression for the total energy of the system: 



N 



(5) 



As it is customary in the analysis of supratransmission, it is convenient to con- 
sider semi-infinite, discrete arrays of nodes, in order to determine the critical thresh- 
old at which the system under study starts to propagate wave signals in the form 
of localized modes. To do that, one may eventually need to consider infinite se- 
quences {un)'^=i of real functions satisfying the differential equations of problem 
^ , for which the sequences of their corresponding derivatives are members of £^ (M) . 
More precisely, one needs to consider sequences (u„)^2 ^ov which J2^=i{'^n)^ is 
convergent at any time. This type of sequences will be called square- summahle. 

Proposition 1 (Macfas-Diaz and Puri [17]). Let {un)n=i be a solution of the system 
of differential equations in f^)- Then, the instantaneous rate of change of the total 
energy of the system with respect to time is given by 



c^(uo - ui)uo - P 



N 



N 



Moreover, if {un)'^^i is a square- summable solution of an infinite class of coupled 
differential equations described by ^3^, then the instantaneous rate of change of total 
energy of the system is obtained by taking the limit when N tends to infinity in the 
right-hand side of 0). □ 

Corollary 2 (Geniet and Leon j21j). Let {un)n=i ^2 a solution of the undamped 
version of f^i- Then, the total energy of the system is given by 



(7) 



E{t) = -c^ f 
Jo 



uo{s){ui{s) - uo(s)) ds. 



□ 



3. Computational techniques 

3.1. A finite-difference scheme. Consider a class of N differential equations 
satisfying (|3]), and let us take a regular partition = to < ti < ■ ■ ■ < tu = T of 
a time interval [0, T], with time step equal to At. For each fc = 0, 1, . . . , Af, let us 
represent by the approximate solution to our problem on the nth lattice site at 
time tk, and the actual value of at the fcth time step by 4>k- Let us convey now 
that 



(8) 





- u^+^ 
















~5ul 
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n — 1 n + 1 

Figure 1 . Forward-difference stencil for the approximation to the 
nth differential equation in ([3]) at time tk, using scheme The 
black circles represent known approximations to the actual solu- 
tions at times tk-i and tk in scheme ([9]), and the crosses denote 
the unknown approximations at time tk+i- 



Then, problem ^ takes the discrete form 
(9) 



s. t. 



m 

'2 



2At 

2 



Slut + TTT^StU. 



■Sut 



0, 



'■N 



7 
2Af 

n't' 

0: 



fc-l 



= J, 



1 < n < A'', 
1< /c < M, 



l<n< N, 
l<k<M, 
l<k< M. 



The forward-difference stencil for this method is presented in Fig. [T] 

In the unbounded situation (that is, when we consider a semi-infinite system of 
coupled Josephson junctions), we choose a large system of N Josephson junctions 
following ([3]), in which the external damping coefficient includes both the effect of 
external damping and a simulation of an absorbing boundary slowly increasing in 
magnitude on the last N — Nq oscillators. Explicitly, we let 7 be the sum of the 
actual value of the external damping coefficient and the function 



(10) 



7'(n) = 0.5 



1 -I- tanh 



/2n-No + N 
\ 6 



where A^o = 50 and N > 200 for our computations. 

It is important to posses discrete schemes to approximate the local energy den- 
sity and the total energy of the problem under study. We associate the following 
expression to finite-difference scheme (O , in order to approximate the value of the 
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n — 1 n + 1 

Figure 2 . Forward-difference stencil for the approximation to the 
nth differential equation in ^ at time tk, using scheme (jl3l) . The 
black circles represent known approximations to the actual solu- 
tions at times tk-i and tk in scheme ([13]), and the crosses denote 
the unknown approximations at time tk+i- 



Hamiltonian of the nth site at the fcth time: 



(11) 



1 

2 

+ ( 







At 










+ l-)2 _ 


2 


2 



,k\2 



fe+1 



Moreover, the total energy of the system is estimated by 



N 



(12) 



/ k 



fe+1 



fe+1 



) + K 



3.2. A second finite-difference scheme. Let = u'^i^^ + 2u^^ -t- u^^^^ . A 

second, finite-difference scheme to approximate solutions to problem ^ is presented 
next. Its forward-difference stencil is depicted in Fig. [21 
(13) 



[Aty 



2 At' 



7 fe 

°-o. 



r2 fe 
,fe+l 



2At 



V{u';+^)-V{ut') 

yk+l _ k-1 



- J = 0, 



1 < n < A, 
1< fc < M, 



u 



1 < n < A, 



s. t. 



cHiu^ - u^) + pStiu^ - u\)IAt = 24>k, l<k<M, 



''N+l 



= 0, 



l<k< M. 



It is worth mentioning that finite-difference schemes ([9|) and ([TS]) are modified 
versions of the one presented in [17j . which in turn is a spatially discretized and 
adapted version of a computational method to approximate radially symmetric 
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solutions of initial-value problems involving a general class of continuous, three- 
dimensional Klein- Gordon equations I16| . 

For approximations obtained via (jl3l) . the Hamiltonian of the nth site of the 
system will be computed using the expression 



,k+l 



(14) 



C 

T 



,fc+l 

''n+l 



,k+l 



J- 



uJi+i -I- ut 



2 2 2 2 

With this notation at hand, the total energy of the system will be approximated 



via 



(15) 



N 



It is important to mention that there are several computational reasons to prefer 
an implicit scheme over an explicit one, in order to approximate solutions to ([3]). 
First of all, it has been noted that some explicit schemes used to approximate 
solutions to modified versions of our problem have proved to be highly unstable 
[12| : on the other hand, the use of an implicit scheme seems inevitable in order to 
approximate the term DfA^Un with a consistency of order at least 0{At)^. 

Also, we must mention that the proposed numerical methods are consistent 
with respect to the problem under study, nonlinear, and require an application of 
Newton's method for systems of nonlinear equations, together with an application 
of Grout's reduction technique for tridiagonal systems pS] , 

Finally, it is worthwhile mentioning here that the stopping criterion used for 
Newton's method in our computations was 



(16) 



max{|u: 



k+i 



1 < n < TV} < 1 X 10" 



4. Numerical properties 

In this section, we prove that finite-difference schemes © and posses prop- 
erties of consistency that make them useful methods in the study of the process of 
nonlinear supratransmission in chains of coupled Josephson junctions. With this 
aim in mind, we present here the most important results concerning the energy 
properties of our methods. 

Proposition 3. The discrete rate of change of energy of a sequence (wj^) satisfying 
finite- difference scheme is given by 
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Proof. We provide a sketch of the proof in view that the task is tedious but, after 
all, a mere algebraic work. Term by term, it is straight-forward to check that 



N 



(18) E 

n 
N 

(19E 



,k+l 



At 



,fc-i 



At 



N 



2 2 



TV 



n=l 



For the sake of convenience, we let 
(21) 

2 ^ 



E 

n=l 

TV 

E 

n=l 

AT 

E 

n=l 



^n+1 



,k+l 



.fc-1 



+K-i-"f;)' 



Using finite-difference scheme ^ and the discrete, boundary conditions imposed 
on the problem, it is straight-forward to check that 

2 ^ 

Bk-Bk-1 = ^ E (MM) ('5*"^;) 
(22) * n_=l 

Notice that the last four terms under the summation symbol form two telescoping 
series, which may be readily simplified. Next, the difference Ek — Ek-i is computed. 
After some simplification, we obtain that 
(23) 



Eh — El 



k-l 



At 



n=l 



2At 



2At 



n—1 ^ ' 



5uq - 5ui\ Stu^ 



2At' 



The result follows now after simplifying the term multiplied by /3, followed by an 
application of the discrete version of Green's First Identity [T7] to the sequence 

Proposition 4. The discrete rate of change of energy of a sequence (w^) satisfying 
finite- difference scheme is given by 



Ek — E, 



N 



(24) 



At 



N 

-7E 



2At 



2At 



')tu1 - dtUQ 
2At 



2Ai 



2At 



Proof. Here, we employ some of the identities used in the proof of Proposition [3l 
Notice first of all that if we define Ck by the formula 
(25) 
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then 
(26) 




We mimic here the proof of Proposition [31 As a consequence of the fact that 
the sequence (u^) satisfies and after an appUcation of the discrete version of 
Green's First Identity [17], it is readily checked that the discrete rate of change of 
the energy between times k and fc — 1 is given by 




whence the resuh foUows after an easy simphfication. □ 

The methods described by finite-difference schemes dH) and ([T3l) are clearly con- 
sistent of order 0{At)^. Moreover, in the following, we give some necessary and 
sufhcient conditions in order for these methods to be stable. Notice that when m 
is a real number, problem (3) represents a discrete, nonlinear modification of the 
classical linear Klein-Gordon equation. 

Proposition 5. Consider finite- difference scheme US]) with V^' = 0, J = and 
a nonnegative constant. Then, At < ^ is a necessary condition for method US]) to 
be linearly stable. 





Proof. We apply the von Neumann stability analysis to method (|13p and we obtain 
the following augmentation matrix: 

(28) A{0 - 
where 

fiO = 2-2c2(Ai)2sin2 

(29) m) = 1 + [c'iAtf + 2/?At] sin^ (0 + ^"'^^'^'^ + and 
MO - l+P(At)^-2/3At]sin^(i) + "-^(^^);"^^^ 

Since all norms in a finite-dimensional space are equivalent and the spectral radius 
p{A) of a square matrix A is the greatest lower bound for the set of matrix-induces 
norms of A, then a necessary condition for linear stability is that p{A) < 1. The 
eigenvalues of A(^) are: 



(30) x^iO- . 
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LU 5 

o . 




2 2.5 3 3.5 4 4.5 5 5.5 

Amplitude 

Figure 3. Graph of total energy versus driving amplitude for a 
system ([3]) subject to a harmonic frequency of 0.8, with coupling 
coefficient equal to 5, and a time period of 10000. Different values 
of external damping were chosen: 7 = (solid), 0.1 (dashed), 0.2 
(dash-dotted), and 0.3 (dotted). 

In order to show that the spectral radius of A{^) is less than or equal to 1, we will 
prove that A±(^) < 1 for all ^ G [0,7r] by considering two cases: 

(1) If A±(0 € C then - mOHO ^ or, equivalently, that 4g(C)/i(0 > 

whence it follows that /i(^) is nonnegative. Moreover, 



(31) 



|A±(e)i 



'MO 



< 1. 



mo V 9(0 

(2) If A±(^) e M and At < i, then we observe first of all that /(^) is non- 
negative with f'^{0 < 4. This case leads to the chain of inequalities 
1 > m > 2 - m > fiO - m- Clearly, 2.g(e) > /(O, so that 



(32) \x±m < 



no + ^Pio - mm) no + \ mo ~ no) 



m) 



< 



m) 



= 1. 



□ 



Proposition 6. Consider finite- difference scheme I113\) with V' — 0, J — and 

a nonnegative constant. This scheme is linearly stable in the infinite norm if 
^ < At < ^. 

7 c 

Proof. A sufficient condition for scheme (|13p to be linearly stable is that 1 1 ^ (C) 1 1 ^ 1 
for all ^ € [0, tt], for some norm || • || ; in our case, we will consider the infinite norm 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Frequency 

Figure 4. Diagram of smallest driving amplitude at which supra- 
transmission occurs versus driving frequency for a system ([3]) sub- 
ject to a harmonic frequency with coupling coefficient equal to 5. 
Different values of external damping were chosen: 7 = (solid), 
0.1 (dashed), 0.2 (dash-dotted), and 0.3 (dotted). 



II • 1 1 00 applied to the augmentation matrix (|28p. Using the same notation as in the 
previous result, observe that the inequality ||v4(^)||oo < 1 is satisfied in case that 



(33) 



9i0 



+ 



HO 



m 



< 1. 



In turn, condition ^ holds if 2 < g{^) + h{^) and 2 < g{^) - h{^). But the former 
inequality is always true, while the latter is satisfied if and only if 



(34) 



1 < 2^Aisin2 ( I 



□ 



Corollary 7. The finite- difference scheme (0) is linearly stable if and only if con- 
dition - < At < - holds. 

7 — — c 

Proof. We apply again the von Neumann stability analysis for the scheme ([9]). The 
augmentation matrix for this scheme is the same as ([25]) when we set /(^) = 2 for 
all 1^ £ [0, tt] and the proof of Propositions [5] and [6] are valid for this case. □ 



5. Application 



Before we provide an application of the numerical techniques presented in this 
work, we must declare that the results of this section were obtained by means of dU 
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Frequency 



Figure 5. Diagram of smallest driving amplitude at which supra- 
transmission occurs versus driving frequency for a system ([3]) sub- 
ject to a harmonic frequency with coupling coefficient equal to 
5. Different values of the Josephson current were chosen: ,7 = 
(solid), 0.1 (dashed), 0.2 (dash-dotted), and 0.3 (dotted). 

and (lisp . Simulations in both the domain of solutions and the energy domain have 
shown that the methods agree excellently. Moreover, it is important to mention 
that we have performed several numerical simulations (not shown here) with both 
methods, for several nonnegative values of /3, 7 and m^, several values of the positive 
constant (? and both V and J equal to zero, observing stability of the methods 
when the corresponding stability conditions are satisfied. 

The validity of the implementation of our methods is checked next against the 
prediction of the nonlinear supratransmission threshold for a semi-infinite, linear 
array of Josephson junctions in parallel and coupled through superconducting wires, 
for which the mathematical model is given by Throughout, we employ a 

coupling coefficient equal to 5, and follow the systematic procedure used in [17) . 
namely: 

(1) For a fixed frequency in the forbidden band-gap, we examine the behav- 
ior of solutions of our problem for different driving amplitudes, using the 
finite-difference schemes proposed in this article. According to the theory, 
there exists a critical amplitude above which the qualitative nature of the 
solutions drastically changes. 

(2) This drastic change in the behavior of the solutions is checked next in the 
energy domain, by making use of the energy schemes associated to each 
finite-difference scheme. A sudden increase in the total energy injected into 
the medium by the driving boundary must be observed above the predicted 
supratransmission threshold. 
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Figure 6. Diagram of smallest driving amplitude at which supra- 
transmission occurs versus driving frequency for a system ([3]) sub- 
ject to a harmonic frequency with coupling coefficient equal to 5. 
Different values of internal damping were chosen: /3 — (solid), 
0.1 (dashed), 0.3 (dash-dotted), and 0.6 (dotted). 

(3) Finally, this procedure is applied to a regular partition {rii}'=o fr^" 
quency interval [0, y/xa^ + 1], obtaining thus a second array {A^j'^Q con- 
taining, for each driving frequency fli , the corresponding critical amplitude 
Ai at which supratransmission starts. In this way, a graph of occurrence of 
critical amplitude versus driving frequency may be constructed. 

(4) This procedure may be repeated for different choices of the parameters of 
the model studied. 

Let us fix a driving frequency 51 = 0.8 in the forbidden band-gap of system ([3]), 
let c = 5, and set all the other scalar parameters of the model equal to zero. Select 
a fixed period of time of 10000, during which the system will be perturbed by the 
harmonic driving (j)(t) = Asin{nt)^ where the value of A will be chosen inside an 
interval that contains the predicted threshold. According to [5], the critical value 
of the continuous-limit case is provided by the expression 



which yields a critical amplitude equal to 3.6 in our case. In these circumstances. 
Fig. [3] presents the graph of the total energy of the system for values of the driving 
amplitude between 2 and 5.5, in the form of a solid line. A drastic increase in the 
total energy is observed around an amplitude of 3.75, in good agreement with the 
continuous-limit case. 

Let us denote by Ei = Ei (A) the total energy of the system described in the pre- 
vious paragraph associated with the amplitude A e [2, 5.5] at time 10000, obtained 



(35) 



As = 2c(l - n^), 
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Figure 7. Diagram of smallest driving amplitude at which supra- 
transmission occurs versus driving frequency for a system ([3]) sub- 
ject to a harmonic frequency with coupling coefficient equal to 5. 
For every € = 0, ±1,±2,±4, the plot Pi corresponds to a system 
with mass m satisfying i/m^ + 1 = 1 + ^ . 



using finite-difference scheme dH). Similarly, let E2 = E2{A) be the corresponding 
energy obtained using scheme (|13p . It is important to mention that our computa- 
tions show that 

(36) msi^{\Ei{A) - E2{A)\ : A £ [2,5.5]} <3x 10"^°. 
Moreover, for amplitudes below the critical threshold, 

(37) msi^{\Ei{A) - E2{A)\ : A £ [2,3.5]} < 8 x 10"^^ 

Next we wish to check the effect of external damping on the occurrence of the 
supratransmission value. To this effect, we fix three different values of external 
damping: 7 = 0.1, 0.2 and 0.3; the results are displayed in Fig. |3l in which 
the corresponding graphs are presented as dashed, dashed-dotted and dotted lines, 
respectively. It is worth noticing that the presence of external damping produces a 
right shift in the occurrence of the supratransmission value and, at the same time, 
a decrease in the total energy of the system. Similar remarks have been made for 
the same problem with Dirichlet boundary data [17] . 

In a next sage, we wish to compute diagrams for the occurrence of supratransmis- 
sion in the presence of all parameters. So, we compute graphs of driving amplitude 
at which supratransmission first starts versus driving frequency for system ^ . Un- 
der the same computational setting as before. Figs. |4] and [5] present such diagrams 
for different values of the external damping coefficient and the generalized Joseph- 
son current, respectively. A quick comparison with similar results obtained using a 
different finite-difference scheme [25] leads us to certify the validity of our method. 
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Additionally, we provide diagrams which evidence the physical implications of in- 
cluding a relativistic mass and internal damping (see Figs. [7] for mass, and [6] for 
internal damping), the results being in qualitative agreement with the Dirichlet 
case, in the following senses: 

• The phenomenon of harmonic phonon quenching still appears in the pres- 
ence of external and internal damping. 

• The discrepancy region due to phonon quenching is shortened as the exter- 
nal damping coefficient is increased. 

• The threshold value at which supratransmission first occurs, for fixed fre- 
quencies outside the discrepancy region, increases as the damping coefficient 
(internal or external) increases. 

• The diagram, for any fixed value of /3 > 0, is approximately equal to the 
corresponding diagram for the undamped system, shifted /3 horizontal units 
to the right. 

• A horizontal shift of vT+rr? — 1 units in the diagram of a sine-Gordon 
system of mass m with respect to the corresponding massless system, is 
observed for small masses and frequencies outside the discrepancy region. 



6. Conclusions 

In this article, we have presented two numerical methods to approximate solu- 
tions of a system of differential equations, which models discrete arrays consisting 
of Josephson junctions coupled with superconducting wires, in which each site is de- 
scribed by a modified, spatially discrete version of the one-dimensional sine-Gordon 
equation. Both methods are consistent and conditionally stable. 

As an important part of the methods introduced, the work provides consistent, 
discrete schemes for the local energies of the links and the total energy of the 
system. These schemes are such that the discrete rates of change of total energy of 
the system are consistent approximations of the corresponding continuous rates of 
change, so that these techniques are ideal methods in the study of the process of 
nonlinear supratransmission. 

As applications, we obtained diagrams of driving amplitude versus driving fre- 
quency, in order to establish the critical amplitude at which supratransmission 
starts. This study was carried out in the presence of the parameters under study, 
that is, internal and external damping, mass and generalized Josephson current. 
Our qualitative results are essentially similar to those obtained for the correspond- 
ing analysis of the Dirichlet case, corroborating thus the validity of our technique. 

Of course, several directions of further research still remain open. For instance, 
the design of computational techniques that employ symplectic methods in the 
analysis of the transmission of energy in discrete arrays of Josephson junctions, is 
a topic that merits close attention. Particularly, one may use numerical methods 
that split the vector field y' as the sum of a conservative contribution fiif) — solved 
in the interval + /i/2] using a symplectic method for the sine-Gordon equation 
[27l|28l [29] — , and a linear dissipative part g{jj) — solved exactly in [i, t + h] with 
an exponential integrator, for instance). The symplectic dynamics may be solved 
then for the next half step, obtaining thus a second-order, splitting method that is 
linearly implicit, and symplectic in the unperturbed case. 
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